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As is widely recognized in Lyapunov analysis, linearized Hamilton's equations of motion have two 
marginal directions for which the Lyapunov exponents vanish. Those directions are the tangent one 
to a Hamiltonian flow and the gradient one of the Hamiltonian function. To separate out these two 
directions and to apply Lyapunov analysis effectively in directions for which Lyapunov exponents 
are not trivial, a geometric method is proposed for natural Hamiltonian systems, in particular. In 
this geometric method, Hamiltonian flows of a natural Hamiltonian system are regarded as geodesic 
flows on the cotangent bundle of a Riemannian manifold with a suitable metric. Stability/instability 
of the geodesic flows is then analyzed by linearized equations of motion which are related to the 
Jacobi equations on the Riemannian manifold. On some geometric setting on the cotangent bundle, 
it is shown that along a geodesic flow in question, there exist Lyapunov vectors such that two of 
them are in the two marginal directions and the others orthogonal to the marginal directions. It is 
also pointed out that Lyapunov vectors with such properties can not be obtained in general by the 
usual method which uses linearized Hamilton's equations of motion. Furthermore, it is observed 
from numerical calculation for a model system that Lyapunov exponents calculated in both methods, 
geometric and usual, coincide with each other, independently of the choice of the methods. 



I. INTRODUCTION 



Natural Hamiltonian systems with many degrees of freedom have Hamiltonian functions of the form 

1 ^ 

In spite of the simple appearance, those Hamiltonian functions having appropriately chosen potential functions are 
used in a wide variety of physical sciences such as plasma physics, condensed matter physics, and celestial mechanics. 
However, the potential functions describe nonlinear interactions, in general, so that chaotic or highly unstable tra- 
jectories take place in respective phase spaces, as is widely recognized. The exponential instability of trajectories are 
measured in terms of Lyapunov exponents, which describe time-averaged properties of chaotic trajectories. Further, 
in the study of directional deviations of chaotic trajectories, Lyapunov vectors will be of great use. 

The Lyapunov exponents and the Lyapunov vectors are defined through linearized Hamilton's equations of motion. 
For the Hamiltonian J^l, the linearized equations take the form 

where X — (Q^, ■ ■ • , , Pi, ■ ■ ■ , Pn) is a 2A^-dimensional vector representing a deviation from a reference trajectory 
{q{t),p{t)) to a nearby trajectory. The linearized equations ||2Jl have 2N linearly independent solutions, which we 
denote by Xa{t), a = 1, • • • , 2N. The Lyapunov vectors Va{t) arc then obtained by orthogonalizing these solutions 
on the Gram-Schmidt method: 

y„.(t) ^ x^it) - 2 7^#PS^^W' « = 1, • • ■ , 2A^, (3) 

^ (Vh(t), Vb(t)) 

where {X, V) denotes the inner product of X and V. The a-th Lyapunov exponent Xa is calculated as 

A.= .»i..,M (4) 

t^oc t \\Va(0)\\ 
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It is to be noted that the values of the Lyapunov exponents are known to be independent of the choice of initial 
values of the Lyapunov vectors except for vanishing Lebesgue measure P, 0] , and that the exponents are ordered as 



Since the Lyapunov exponents are time- averaged quantities, they are suitable for the study of statistic properties 
of Hamiltonian systems. For example, phase transitions are investigated by the use of Lyapunov exponents. In fact, 
the second-order phase transition 3] and the Kosterlitz-Thoulcss transition J| are characterized by the discontinuity 
in the largest Lyapunov exponents and by a sudden change in the gradient of the largest Lyapunov exponent against 
energy, respectively. Further, the sum of all positive Lyapunov exponents, which is also viewed as a function of energy, 
is used in the discussion of a dynamical phase transition 0, according to which trajectory's phase transition from 
nearly-integrable behavior to chaotic behavior occurs in an energy region in which the sum of positive exponents 
breaks into a rapid increase against energy. In contrast with this, the Lyapunov vectors are expected to be useful in 
studying dynamical behavior of chaotic trajectories, since they serve as time series 7]. 

Suppose that a reference trajectory is given in a phase space. Then, according to ((SJ, one can form 2N linearly 
independent Lyapunov vectors from solutions to the linearized equations of motion along the reference trajectory. 
However, two of the Lyapunov vectors which are associated with Lyapunov exponents, Ajv and Ajv+i, are considered 
as marginal, since Aat and Atv+i should vanish, as is widely recognized. One of those two Lyapunov vectors is 
the tangent vector to the trajectory, Xh, and the other the gradient vector of Hamiltonian function, gradiJ. We 
may interpret these vectors as follows: The displacement in the direction X.h is regarded just as a certain time 
displacement in the reference trajectory, and the displacement in the direction gradif will give rise to a transfer to 
a nearby trajectory with an energy value different from that of the reference trajectory. In view of this, in order to 
analyze the instability of trajectories, we are allowed to require that the two directions pointed by the vectors Xu and 
gradi? be separated out from the other 2N — 2 directions. Put another way, the requirement means that a Lyapunov 
vector which is orthogonal to the plane spanned by Xh and gradiJ at an initial instant has to be orthogonal to 
the plane spanned by Xh and gradi/ at every instant. If the requirement is fulfilled, we will be able to discuss the 
instability of trajectories without influence of the two marginal directions. 

Unfortunately, the usual method of Lyapunov analysis on the basis of the equations Q does not satisfy the 
requirement in general. This is because for any solution X{t) to Q one has 



so that one obtains {X^gv&AH) = at any instant if (X,gradiJ} |f=o = at an initial instant, but, in general, by no 
means one can make {X,Xh) vanish at any instant, so that even the first Lyapunov vector, Vi, can not be made 
orthogonal to the plane spanned by Xu and gradiJ at every instant. 

A way to construct Lyapunov vectors which satisfy the above-stated requirement is to adopt linearized equations 
of different type from the usual one (0). To take a geometric approach to Hamilton's equations of motion is a step 
toward finding such Lyapunov vectors. As for the geometric approach, it is known that if the total energy of the 
natural dynamical system is fixed at E, Newton's equations of motion can be equivalently expressed as geodesic 
equations on a Riemannian manifold (M, gij ) , where M is a subspace of the configuration space defined by 
M = {q£ R^\E - V{q) > 0} and gij is the Jacobi metric defined by gij{q) — 2(E — V{q))Sij. Then the linearized 
equations of the geodesic equations are given by the Jacobi equations of the form 



Ai > A2 > • • ■ > A2Ar. 



^(X,gradH) =0, 



(5) 




N 



(6) 



where Rji^/ are the components of the Riemann curvature tensor, and s is the arc length defined as 




(7) 



On the other hand, the linearized Newton's equations are put in the form 




N 



(8) 



Equations ^ and are not transformed to each other through the parameter transformation {7|), while Newton's 
equations of motion and geodesic equations for the Jacobi metric are transformed to each other. This geometric 



method has been introduced in the estimation of the largest Lyapunov exponent Ai with the aid of statistical mechanics 
1^ lE, il^ ilijj • They studied instability of geodesies through the Jacobi equation, a second-order differential equation, 
while the Lyapunov analysis needs first-order differential equations. 

The geometric approach we will take in this article is to be made on the cotangent bundle T* M of the Riemannian 
manifold M in order to find first-order differential equations associated with © and thereby to construct Lyapunov 
vectors which satisfy the above-stated requirements. We will first work with generic linearized Hamilton's equations 
of motion on T*A/, and then specialize the resultant equations to linearized Hamilton's equations for geodesic flows 
on T*M , which will be found to project to the Jacobi equations on M . Further, we will introduce a lifted metric on 
the cotangent bundle T*M to make it possible to discuss the orthogonality of vector fields on T*M. The lifted metric 
may be called the Sasaki metric. On this setting, we will be able to find Lyapunov vectors satisfying the above-stated 
requirements along any geodesic flow on T*M. Put in detail, it will be shown that along any geodesic flow on T*M, 
there exist Lyapunov vectors such that those associated with the vanishing Lyapunov exponents Xn and Aat+i are 
Xh and gradJf, respectively, and the other 2N — 2 Lyapunov vectors arc all orthogonal to the plane spanned by Xh 
and gradT? at each point of the geodesic flow. 

This article is organized as follows: Section ITll cont ains a brief review of geodesies and Jacobi fields and, in particular, 
of the Jacobi metric, whose geodesies are equivalent to trajectories of the natural dynamical system with a fixed total 
energy. In Sec^l and succeeding sections, Einstein's summation convention is adopted, and we choose to denote 
by (x*) local coordinates on a general m-dimensional Riemannian manifold, and by (q^) the Cartesian coordinates 
on . Section Hill is concerned with geodesic flows on the cotangent bundle T*M, which project to geodesies on 
M . To describe geodesic flows in a more geometric way, we introduce an adapted frame and a lifted Riemannian 
metric on T*M. Linearized Hamilton's equations of motion are discussed in Section Hvl and it will be shown that 
there exist Lyapunov vectors which satisfy the above-stated requirement along a geodesic flow on T*M. Section Ivl 
is for numerical calculations for a model system with three degrees of freedom. Lyapunov vectors and Lyapunov 
exponents are calculated numerically for the model system in both geometric and usual methods to compare the 
respective results. It will be shown that Lyapunov exponents calculated in respective methods coincide with each 
other, independently of the choice of methods. Section lVll is devoted to concluding remarks. Appendices are attached 
in which related topics on geometry of cotangent bundles and a symplectic implicit Runge-Kutta method for numerical 
integration are reviewed. In particular, lifting vector fields on M to T*M and the Levi-Civita connection with respect 
to the lifted metric on T*M are discussed. 



Let (M, g) be an m-dimensional Riemannian manifold with metric g. The metric induces the Levi-Civita connection 
V on M; for vector fields Y, Z ^ SB{M), SG{M) denoting the set of vector fields on M, the covariant derivative 
WyZ is defined, in terms of local coordinates {x^, • • • , x™), to be 



II. GEODESICS AND JACOBI FIELDS 



A. Jacobi equations 




where (y ) and (Z*) are components of Y and Z, respectively, the Christoffel symbols Tjj, are defined as 




with components of the metric 





For a geodesic c(s) with s the arc length parameter, the tangent vector ^ to the geodesic satisfies the geodesic 
equation 



4 



where 



ds dx^ 



c(s) 



with 

ds^ = g.jdx'Ax^. (10) 

We are interested in stabiUty/instabihty of geodesies. To this end, we consider a congruence of geodesies which 
looks hke a fluid whose flow lines are geodesies with the c(s) as a member of them. Then we may consider that the 
tangent vector ^ to c(s) is extended to be a vector field defined in a neighborhood of the original geodesic c(s). We 
may also assume that there exists a vector field Y satisfying the condition 

[tY]^Q (11) 

in the same domain as that for ^. The condition (|ll|l means that a geodesic with ^ its tangent vector is carried 
congruently to another infinitesimally nearby geodesic by the infinitesimal transformation Y . Thus Y is viewed as a 
deviation of geodesies. The vector field ^ may have singularity at which ^ is not defined uniquely, and Y may vanish 
there. With this in mind, we operate V^^ = with Vy and use the definition of the Riemann curvature tensor and 
of the torsion, which vanishes identically, to obtain the Jacobi equation 

V^y^Y + R{Y,m = 0. (12) 

Here, as is well known, the torsion tensor and the Riemann curvature tensor are defined, respectively, to be 

T{Y., Z) = VyZ - VzY - [Y, Z], 

R{Y, Z)W = Vy^zW - Vz^yW - V[y^z]W, 

where Y , Z,W € cB{M), and the Riemann curvature tensor has symmetries such as 

9iR{Y, Z)W, U) = -g{R{Z, Y)W, U) 

(13) 

= -g{R{Y, Z)U, W) = g{R{W, U)Y , Z), 
where Y , Z ,W ,U e S6'{M). Local components of the Riemann curvature tensor are expressed as 

R -H -n(n( ^ ^ \ ^ ^ 



dx"^ ' dx^ ) dx^ ' dx^ 

r> I _ jk _ <J^ ik -pi pm p£ pm 
dx^ dx^ ^ im^ jk ^ jm^ ik- 

In the next subsection, we will give an example of Riemannian metrics whose geodesies are equivalent to trajectories 
of Newton's equations of motion 

d\' dV , , 

df+a?=0' -V-,A^. (14) 

Then, in order to analyze stability/instability of trajectories of the natural Hamiltonian system, we can deal with the 
Jacobi equation, a linearization of the geodesic equation. However, the Jacobi equations in their original form are not 
suitable to Lyapunov analysis. 

B. Geodesies for the Jacobi metric 

Consider equations of motion, Ea. (|14|l . on R^ , which we call a natural dynamical system. Let Mj be an open 
submanifold of the configuration space R^ , which is defined to be 

Mj = {qeR''\V{q)<E}. (15) 
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As is well known, if energy is fixed at E, almost all trajectories are confined in Mj when N > 2. On the other hand, 
the Jacobi metric g j is defined, in Mj, to be 



(gj),, = 2[E - y(g)]<5. 



(16) 



According to Maupertuis's Principle of Least Action, an extremal of the action, the integral of the kinetic energy 
along possible paths, provides an actual trajectory of total energy E. This principle can also be stated as follows: An 
extremal of the variational problem of lengths of paths with respect to the Jacobi metric provides an actual trajectory 
of the total energy 

HmFl- From Eq.^ along with the Jacob i metric g j, the arc length parameter s is shown to 
be related to the time parameter t by 

ds^ = 4[i;-y((?)]2dt2, 
and the tangent vector to a geodesic is always unity accordingly: 

dq* dq^ 



[2iE~V)Y 



(17) 



Since the Christoffel symbols for the Jacobi metric (|16|l are given by 

-1 



2{E - V) 



dV 



the geodesic equations for the Jacobi metric are expressed as 

1 dV Aq^ dg' 



dV 
ds2 



dV 



E-Vdqi ds ds A{E-VYdq' 



0, 



(18) 



which prove to be equivalent to Newton's equations of motion (|14|l on account of (|17|l . However, the Jacobi equations 
with the curvature tensor for the Jacobi metric are not brought into the same equations as ((SJ, a linearization 
of Newton's equations of motion, in general. Components of the curvature tensor for g j are indeed put in the form 



Rijki = CuSjk + CjkSii — Cik6j£ — CjiSik, 



where 



dV dV 



1 



~ dq^dq^ ^ 2{E -V)d^W ^ <E - V) d^W'^' 



III. GEODESIC FLOWS ON COTANGENT BUNDLES 



In the previous section, we have mentioned that trajectories of a natural dynamical system with a fixed energy 
can be regarded as geodesies on a suitable Riemannian manifold, and that stability/instability of the trajectories are 
analyzed through the Jacobi equation, a linearization of the geodesic equation. However, the Jacobi equation is a 
second-order differential equation, while Lyapunov analysis is applied to first-order differential equations. We hence 
need a first-order differential equation associated with the Jacobi equation in order to apply Lyapunov analysis. To 
find such a first-order differential equation, we are working on the cotangent bundle T*M of a Riemannian manifold 
M along with some geometric setting-ups on T*M. Related topics on T*AI will be described in Appendix IXI 

At first, let us be reminded of a minimum on cotangent bundles. Let M be an m-dimensional Riemannian manifold 
endowed with Riemannian metric g = gijdx^ O dx^ , and T*M the cotangent bundle of M with the projection 
TT : T*M M. Let (cc^) and (x*,pi) be local coordinates in an open subset U C M and in tt~^{U), respectively. 
Further, let (x',|5j) be another local coordinates in Tr^^{U) with tt^^{U) D Tr^^{U) ^ 0. Then one has the coordinate 
transformation in the intersection tt~^(U) fl tt~^{U), 



_■ _ dx' 



(19) 



6 



A. Geodesic flows 



We recall that Newton's equations of motion have been already "geometrized" so as to be geodesic equations on 
a suitable Riemannian manifold, so that further external force doesn't need to be taken into account anymore. In 
other words, we have only to consider a free particle motion on M . In the Hamiltonian formalism, the Hamiltonian 
we then have to study should be given, on T*M, by 

K^^g^^{x)p^pj, (20) 

where '.— {gij)~^- Hamilton's equations of motion for K arc then put in the form 

dx' dK 

as dpi ^21) 



dpi dK „ 

d7 = -a^^^^-^'^^^' 



where use has been made of the equality 



QgM 



It is an easy matter to show that Ea. (|21|) projects to geodesic equation on Af. In fact, put together, differentiation 
of the first equation of (|?T|l with respect to s and the second equation of (j?T)) along with the above equality provide 
geodesic equations. As is well known, Ea. (|21|l is associated with the Hamiltonian vector field Xk given by 

dK d dK d d fe, . d 

Integral curves of H22|) are called geodesic flows. We note here that the nomenclature "geodesic flows" are usually 
assigned to the corresponding flows on the tangent bundle, but we use the word for convenience' sake. 

It is worth noting here that how geodesic flows on T*A'I project to geodesies on M. Let P{s) — {x{s),p{s)) be a 
geodesic flow with an initial value (x(0),p(0)) = (a, b) with g^^hthj = 1. Define a tangent vector v to M at a = 7r(P(0)) 
by = Then the projection x{s) = ■n{P{s)) is a geodesic with the initial value (a;(0),i(0)) — [a,v). Varying 

h G T*M with g^^bibj — 1 but fixing a, we obtain an (m — l)-parameter family of geodesic flows on T*M, which 
projects to an (m — l)-parameter family of geodesies passing the point a of M. Furthermore, the vector field Xk 
projects to the tangent vector field ^ to this family of geodesies M. However, ^ has singularity only at o € M in a 
neighborhood of a. We have to note that if all geodesic flows on T*M project to geodesies M, those geodesies may 
not define such a tangent vector field uniquely. If there is another (m— l)-parameter family of geodesic flows on T*M, 
it may project to another (m — l)-parameter family of geodesies on AI along with a tangent vector field like ^. 



B. Adapted frames 

To describe geodesic flows in a more geometric way, we introduce an adapted frame and an adapted coframe on 
7r~^(C/) C T*M by the use of the Christoffel symbols r*^, on M. The adapted frame and coframe are defined, in 
7r^i(f7), to be 



and 



dx\ ^ dp^ - pkT'ljdx^ , (24) 



respectively, where i — i + m. These frames are dual to each other, i.e., satisfy 
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If there is another adapted frame {Z^i, 13-} in an open set TT ^{U) and if the intersection tt ^{U)rnT ^{U) is non-empty, 
then from H19|l it follows that the adapted frames are subject to the transformation 

For adapted coframes, an analogous transformation holds as well. 

The transformation H25|l implies that Di, i = 1, • • ■ , m, and Dj, i = m + 1, • • ■ , 2to, define, respectively, subspaces 
Hp and Vp of the tangent space TpT*M at each point P € T*M independently of the choice of adapted frames. 
Thus one obtains a direct sum decomposition of the tangent space to T*M at each point P G T*AI; 

TpT*M ^ Hp®Vp. (26) 

The subspaces Hp and Vp are called the horizontal and the vertical subspace of TpT*A{, respectively. We notice here 
that Hp and T^^i^p-^M are isomorphic as vector spaces. Note further that the transformation rule for the standard 
frame {d/dx^, d/dpi] is mixed up, so that one cannot define a subspace, say, span{9/9a;'} independently of the choice 
of natural frames. See Q| for adapted frames on the tangent bundle TM. 

In terms of the adapted frame, the Hamiltonian vector field Xk becomes expressed as 

Xk = Dj{K)D, - D,{K)Dj - g'^PjD,, (27) 

which shows that Xk is a horizontal vector field and further that geodesic flows are horizontal curves in the sense 
that the tangent vectors to them are always horizontal. 



C. The Sasaki metric 



As is already seen, the tangent space to T*M at each point of T*M is decomposed into a direct sum. We can define 
a metric on T*M so that the decomposition may be orthogonal direct sum. One of such metrics is the Sasaki metric, 
which is a lifted metric g given by 

g = g,^9'®e' +g'^e'®e^. (28) 

This metric is defined independently of the choice of adapted coframes. We notice here that the Sasaki metric was 
first introduced on the tangent bundle TM but we use the same nomenclature on the cotangent bundle T*M as 
well. 

By using the Sasaki metric, the arc length on T*M is defined as 

da^ = g.jdx'dx^ + g'^idp, - pferf„dx")(dp, - p^F^'^dx"). 

It then turns out that geodesic flows on T*M take the same arc length as the corresponding geodesies on M have, since 
one has dcr^ = gijdx^dx^ = ds^ for horizontal curves, and since geodesic flows are horizontal. Hence the parameter s 
used in Hamilton's equation lf?T|l may be interpreted as the arc length on M, so that the geodesic x{s) = tt{P{s)) on 
M is described in the arc length parameter. 

We will adopt the Sasaki metric on T*M to discuss orthogonality of Lyapunov vectors on T*M in the next section. 



IV. LYAPUNOV ANALYSIS OF GEODESIC FLOWS 



On the basis of the geometric setting-up, we are to find a first-order differential equation associated with the Jacobi 
equation, and thereby discuss Lyapunov vectors. 



A. Linearization of Hamilton's equations of motion 



For a general Hamiltonian function H , linearized Hamilton's equations of motion are put, as is well-known, in the 
form 

di-' _ d'H ^ d'H ^-^ 
ds dpidx^ dpidpj ' 

dX~ _ __d^p _ d^H ^- 
ds dx'^dxi dx^dpj ' 
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where X = X'^di + X'^d^ stands for a deviation of Hamiltonian flows, where di = d/dx^ and 9j = d/dpi. These 
equations can be obtained from the condition [X , Xfj] = as weh, where Xu is the Hamiltonian vector field, 



X 



dH d dH d 



H 



dpi dx^ dx^ dpi 

In fact, the condition [X,Xif] = restricted to a prescribed Hamiltonian flow P{s) — {x{s),p{s)) provides 



[X,Xh]\p 



is) 



X^ 



dpidxi 



X' 



dptdpj 



dx^dxi 
dX' d 



ds dx^ 



X^ 



Pis) 



dx^dpj 
dX~ d 



d 

dx^ 
X^ 



P(s) 

_d_ 

dpi 



P(s) 



ds dp^ 



Pis) 



where we have used the formula 



dX' dH dX' dH dX' 



ds dpk dx^ dx^ dpk 



Xh{X' 



and a similar formula for dX*/ds. It is to be noted here that the condition [JC, JC^f] = implies that a Hamiltonian 
flow, an integral curve of Xh, is dragged to another infinitesimally nearby Hamiltonian flow by the inflnitesimal 
transformation X, i.e., X is a, deviation of Hamiltonian flows. With this in mind, we can obtain linearized equations 
with respect to the adapted frame, if we calculate [X,Xh] = with X and X^ expressed as 

X = X'D, + X'Dj, 
Xh =Dj{H)D,-D,{H)Dj, 

respectively, and restrict the resultant equation to a prescribed flow P{s). We note here that the components (X', X'-) 
with respect to the adapted frame transform according to 



dx^ ox^ 



(29) 



A long but straightforward calculation of [X,X/f]|p(s) = then provides linearized Hamilton's equations of motion 
as follows: 



dX^ 
~d7 

dX' 
~di 



d^H 



d^H 



dpidxi dpidp, 
d^H d^H 



PkTl 



X^ 



d^H 

dpidp-j 
( d^H 



X\ 



d^H 



dx^dx^ dx^dp. 



^^^^'^ + [d^ + a^^"^- ' ^^^^^ 



(30) 



dx^ 9a;™ dpm 



X' 



d^H 



d^H 



dx'^dpj dpfdpj 



XK 



where use has been made of the formula 

dX 



ds 



D-{H)D,{X^) - D,{H)Dj{X^) = Xh{X% 



and of a similar formula for dX^ /ds. 

In what follows, we take the Hamiltonian given by Ea. H2U|l . The equation of deviation ()30|l then takes the form 



dX' 
~d7 
dX~ 
~d7 



-T)^g^'ppX' +g'^X\ 
-RjkH9''''Pn9"'PhX^ + rl.g^'piX 



(31) 
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The right-hand-side of Ea. pifl must be evaluated along a geodesic flow P{s) = {x{s),p{s)). Since one has g'^^pj{s) 
da;* 

— — =: C{s) along the geodesic flow, Eg. (13 111 can be brought into the form 
ds 



(32) 

We can show that this system of equations is the first-order differential equation which project to the Jacobi equation 
and hence may be called the lifted Jacobi equation. The proof runs as follows: On account of (|29|l . the quantities 
(X*(s)) and (X'(s)) may be viewed as a tangent and a cotangent vector to M along the geodesic x{s), so that the 
first equation of 1)32(1 . rewritten as 

d X* - 

implies that {g^^ X^{s)) is equal to the covariant derivative of (X*(s)) along the geodesic x(s). The second equation 
of H32|l then implies that 

Ayi - 

— +T),eY^=-R^,,^e^'X^ with Y^^g^^X^. 

The above two equations are put together to yield the Jacobi equation for Yx — {X^{s)), 
where stands for the covariant derivation along the geodesic x{s). 



B. Lyapunov vectors 

Here we show that solutions to Eg. 1(32(1 satisfy the requirement stated in Introduction in the Haniiltonian system 
with the Hamiltonian K given in ((20(1 . As for the gradient of K, we note that the differential dK and the gradient of 
K, gradi^, are put in the form 

dK = g'^'pkO^ gradX = p,Dj, 
respectively, where the gradient of a function F on T*M, gradF, is defined through 

g(gradF,X) =dF(X) 

for any vector field X e 3&{T*M). 

It is an easy matter to verify that Eg. ((32(1 is satisfied by Xk, the tangent vector to a Hamiltonian flow P{s) or 
a geodesic flow in T*M. In fact, the tangent vector Xk to P(s) is given l(?7|l . and has the components, X^{s) = 
g^^Pi{s), X'^{s) = 0, satisfying 1(32(1 . While the gradient vector along the Hamiltonian flow P(s), which is denoted 
by gradif (s) for simplicity, is not a solution to the linearized equation 1(32(1 . the vector gradi4r(s) -I- sXk{s) — 
Pi{s)Dj + s g'^^pk{s)Di is a solution to 1(32(1. as is easily verified. Taking this into account, we wish to decompose the 
tangent space Tp(^syr*M to T*M at every point P(s) of a geodesic flow into the direct sum of the plane spanned by 
both Xk{s) and gradiir(s) and the subspace transversal to the plane. Let us define subspaces Np(^s-^ and Ep(^s-^ to be 

7Vp(,) = {X e Tp^,^T*M I X = aXK{s)+pgTadK{s), a, /3 £ R}, 

(33) 

£;p(,) - {X e Tp^s)T*M I g{X,XK{s)) = 0, g(X, gradi^(s)) = 0}, 

respectively, where i?_p(s) is the orthogonal complement of A^p(s) with respect to the Sasaki metric g. Thus we have 
the orthogonal direct sum decomposition. 



(34) 
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We wish to show that these subspaces are invariant under any solution to the hnearized equation (|32|l . To this end, 
we have to verify: 

Theorem: A solution X{s) to the linearized equation (|32|l which is in A'^p(o) (resp. in £'p(o)) at an initial moment 
s = keeps belonging to A^p(s) (resp. to Ep^g-^) at any instant s. 

The proof of this statement is carried out as follows: As we have already shown, Xk{s) and gradA'(s) + sXk{s) 
are solutions to H32I) . so that the linear combination of them, a Xk{s) + P{gra.dK{s) + s Xk{s)) = {a + Ps)Xk{s) + 
/3grad-fC(s), is also in A^p(s) at any instant s, which proves the invariance of A^p(s) under the linearized flow X{s). 
To prove the invariance of Ep(^s), we consider the temporal evolution of g{X, Xk) with X a solution to H32|l . We are 
to show that 

^g{X,XK)=dK{X), 

2 (35) 

^g(X,XK)-0. 
ds^ 

We can carry out the proof of these equations in the manner of mechanics as follows: Note that g(X, Xk) = (^{X), 
where 9 is the standard one-form on T*M, i.e., 9 = Pidx^ in local coordinates. Then differentiation of 9{X) with 
respect to s results in 

^9{X) = CxA9iX)) 
as 

= {Cx,MX) + 9{[Xk,X]) 
- idi{XK)9 + iiXK)d9)iX) 
= {di9{XK))-dK){X) 
= dK{X), 

where use has been made of (i) the definition of the Lie derivative of one- forms, (ii) the condition [X, X^] = 0, (iii) 
the Cartan's formula for the Lie derivation, (iv) i{XK)d9 = —dK, and (v) the equality 9{Xk) — 2A' due to the 
homogeneity of K in pi. Thus we obtain the former equation of (|35|l . Differentiating the former equation of (|35|l with 
respect to s using the equation 

^dK{X) ^ ^g{gradK,X) = 0, 
ds ds 

a similar equation to we obtain the latter of (|35|l . Now Ea. (|35|l is integrated to give 

g{X,XK)\p{s) =g{X,XK)\p(o) + sdK{X)\p(^o)- 

Since dK{X) = g{X ,gTa,dK), the above equation implies that X{s) G £'p(s) if X(0) G -Ep(o)- This ends the proof of 
the invariance of Ep(^gj under the linearized flow X{s). ■ 

On the basis of the decomposition (|34|l . we can construct a set of Lyapunov vectors {Va}, a — 1, - ■ ■ , 2m, satisfying 
the requirement mentioned in Scc^ We are thinking of the Riemannian manifold {Mj,gj) introduced in Sec III Bl 
and hence m = N . The first — 1 linearly independent solutions, {Xa(s)}, a=l,---,A''— 1, to the lifted Jacobi 
equation H32(l are chosen in £'p(s), which are orthogonalized to give first A^ — 1 Lyapunov vectors {Vi, • • • , V m-i}- 
The A^-th and (A^-l-l)-th Lyapunov vectors are chosen in Afp(s) so as to be V m{s) = Xx{s) and V m+i{s) — gradAr(s), 
respectively. This is because they are mutual orthogonal and because Xk{s) and gradAr(s) -I- s Xk{s) are solutions 
to the linearized equation and further orthogonal to the first A^ — 1 Lyapunov vectors. Note in addition that any 
solution X[s) staying in Npi^^-^ becomes asymptotically parallel to Xk{s) as s ^ oo, so that Xk(s) is assigned to the 
A^-th Lyapunov vector and gradi?(s) to the (A^-|- l)-th Lyapunov vector, respectively. The remaining A^ — 1 Lyapunov 
vectors are chosen in Ep(^s) , which are orthogonal to Xk and gradAT as well as to the first A^ — 1 Lyapunov vectors 
by the very definition. Consequently, from solutions to H32|l with the initial values chosen so as to satisfy 

(a) Xjv(O) = Xx(0), Xjv+i(0) = grad^C(O), 

(b) X„(0) ± {Xk{s), gradi^(O)}, a = 1, • • • , TV - 1, A^ + 2, • • • , 2A, 



at the initial moment s = 0, we can obtain expectedly a set of Lyapunov vectors such that 
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(i) Vn{s) = Xk{s), Vn+i{s) = gradi^(s), 

(ii) Va{s) ±{XKis),gTeidK{s)}, a=l,--- ,N -l,N + 2,--- ,2N. 

From the property (i), we can observe that the Lyapunov exponents Xn and \n+i vanish indeed. In fact, since 

giXK,XK) = gigradK, gradX) = 2K 
is constant along any geodesic flow, one has Xn — Xn+i — from the formula 

V. NUMERICAL CALCULATIONS FOR COMPARISON 

In this section, we are to compare the geometric method and the usual method through a model system with 3 
degrees of freedom, by numerically calculating Lyapunov exponents and Lyapunov vectors in respective methods. We 
will find that the Lyapunov exponents calculated in respective methods coincide with each other, independently of 
the choice of methods, while the Lyapunov vectors calculated on respective setting-ups exhibit different behaviors to 
each other, depending on the method chosen. 

A. Comparison of setting-ups in respective methods 



TABLE I: Comparison between the usual method and the geometric method. The A'^-dimensional manifold AIj is defined in 
Eg.lT^. and is the Euclidean metric. Note that T*Mj = Mj x R'^ . 

Method Configuration Phase Metric Hamiltonian Linearized 
space space equation 



Usual Mj Mj X H [Eq.((TJ] Eq.0 

Geometric Mj T*Mj Qj K [Eq.llUll] Ea. lf^ 



For a natural Hamiltonian system with N degrees of freedom, setting-ups for Lyapunov analysis both in the 
geometric method and in the usual method are summarized in Tab^ We note here that the metric introduced on 
the phase space Mj x R in the usual method is, of course, the Euclidean metric defined, as usual, to be 

Qe ^ Sijdq^ ® dq^ + S^^dpi ® dp-j. 

As was pointed out in Sec III Bl the geodesic equations for the Jacobi metric are equivalent to Newton's equations 
of motion for a natural dynamical system with energy E. We now verify this fact in the Hamiltonian formalism. The 
Hamiltonian vector fields Xk and Xh which are defined on the same phase space in respective manners are given by 



X..,>,A, (36) 



respectively, where (/'^ = 5^^ /2{E — V). A straightforward calculation along with H18|) and \ ^Pi + V = E then 
provides 

= 2{E^-V)^"' ^^^^ 

which implies that Hamiltonian flows both in the geometric method and in the usual method coincide within the 
change of parameters, ds/dt = 2{E — V{q)). Thus, along the same flow (up to the parameter change), we can 
compare numerically tangent vectors such as solutions to linearized equations of motion and Lyapunov vectors. In 
the following, X^^\s) and X(t) denote solutions to the linearized equations of motion in the geometric method and 
in the usual method, respectively. 
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B. Orthogonal relations in the usual method 

In Sec II VI we have shown that Lyapunov vectors in the geometric method can be chosen so that two of them may 
be the tangent vectors to the Hamihonian flow in question and the gradient vector of the Hamiltonian function along 
the flow, and the others be orthogonal to those two vectors. In this subsection, we remark that such orthogonal 
relations holds for part of Lyapunov vectors even in the usual method, in which the Euclidean metric is adopted 
in Mj X R^. 

Let Xi{t)^ ■ ■ ■ ,X2Ar(t) be linearly independent solutions to Eq.(|2Jl, for which the initial conditions are taken in 
such a manner that 

(a) Xn{0) = Xh{Q), Xn+i{0) = gradi/(0), 

(b) X,(0) ± {XniO), gradi?(0)}, a = 1, • • • , TV - 1, iV + 2, • • • , 2iV, 

where Xh and gradi? are the Hamiltonian vector field for H and the gradient vector field of H, respectively. Let 
Vi{t), ■ ■ ■ , V2N{t) be Lyapunov vectors formed from Xa{t), a — 1, • • ■ , 2N. Then the following two properties hold 
true: 

(i) Vi{t), - ■ ■ , yjv(t) are always orthogonal to gradi7(t), 

(ii) VN+i{t), ■ ■ ■ , V^2w(0 are always orthogonal to X// (t). 

The property (i) is easily shown to hold from Eq.©. In fact, solutions Xi{t), - ■ ■ ,Xf^{t) to the linearized equations 
(O are always orthogonal to gradi?(t), if they are initially orthogonal to gradiJ(O). Hence the Lyapunov vectors 
Vi{t), • • • , VAr(t) are always orthogonal to gradi?(t), since the A^-dimensional space spanned by Vi{t), • • • , VAr(i) is 
the same as that spanned hy Xi(t), ■ ■ ■ , X i^(t). For the proof of the property (ii), we use the fact that the Hamiltonian 
vector field Xnit) is a solution to the linearized Hamilton's equation (O, so that one has XN{t) = Xnit). Then 
Xnit) is in the A^-dimensional space spanned hy Xi(t), ■ ■ ■ , X N{t), and hence in that spanned by Vi{t), ■ ■ ■ ,V N{t). 
By definition, the Lyapunov vectors VN+i{t), ■ ■ ■ , V2N{t) are orthogonal to Vi{t), ■ ■ ■ , VN{t), and hence to Xnit). 

The above two properties will be confirmed as well by numerical calculations for a model system in a later subsection. 
Moreover, by numerical calculations in the usual method, we will observe that V ]\[^2it), • • • , V2Ni't) are not always 
orthogonal to gradi/(t), and that Vi{t), • • • , VN-iit) are not always orthogonal to Xnit), either. We recall here 
that, in the geometric method, Lyapunov vectors V'j^l|_2(s), • • • , ^2^(5) are always orthogonal to gra,dK{s), and that 
v[^\s),--- are always orthogonal to Xk{s), which will be confirmed as well by numerical calculations 

for the model system. Here, by V^^^ and Va, we denote the Lyapunov vectors which are obtained in the geometric 
method and in the usual method, respectively, to tell the difference between them. 

C. Initial conditions 

To compare numerical computation results calculated both in the geometric method and in the usual method, 
we have to set both Hamilton's equations of motion and linearized equations of motion to share the same initial 
conditions. Hence, in particular, we come to require that the initial conditions for linearized equations of motion 
are taken to be subject to the conditions (a) and (b) mentioned in Sec lVBI in the usual method as well as in the 
geometric method. In this subsection, we discuss how one can set such initial conditions, in spite of the difference 
between metrics used. 

We take a number of initial values, P{0) — (9^ (0),pj(0)), for Hamiltonian flows on T*Mj in such a manner that 
grady vanishes at the initial point P(0), where grad^ is defined with respect to both the Euclidean metric and the 
Jacobi metric on the configuration space Mj, but the equation gradV^ = defines the same points, independently 
of the metric chosen. Since the phase spaces in both methods are in common, and since Hamiltonian flows in both 
methods are also in common up to the change of parameters, we will obtain a number of Hamiltonian flows in common 
after integration. We have also to note that the condition gradF = at the initial point implies that the Christoffel 
symbols T^-^'s defined by lfTH|l vanish also there, so that the Jacobi metric is put, at the initial point, in the form 

9,/|p(o) = 2Wo5^,dx' ® Ax' + {2Wo)-H''dp^ ® dpj, (38) 

where Wq = E - V{P{Q)). 

Let X^^'^ and Xa denote solutions to linearized equations H32|) and (|2Jl, respectively. 

According to the procedure stated in Sec lIVI initial conditions for linearized equations in the geometric method 
are set as follows: 
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(a) X^^\0) = Xk(0), x'^^l.iO) = gradif(O), 

(b) Xis)(0)ei?p(o)niJp(o), a=l,-- - ,7V-1, 
Xi''' (0) e Ep(^o) n Vpf^o) , 6 = TV + 2, • • • , 27V. 

See (|26ll and (|33|) for the definitions of -ffp(s), Vpf^), and Epf^gy 

Initial conditions for the hnearized equations IpJ in the usual method are set as 

Xa{0)^X[3\Q), a=l,...,27V. 

We here have to verify that these initial vectors Xa(0), a = 1, • • • , 2N , are indeed subject to the initial conditions 
(a) and (b) stated in Sec lVBI The verification of this is carried out as follows: By definition, one has Xjv(O) — 
Xk{Q), and further Xk{Q) = Xh{Q)I2{E ~ V{q{Q))) from so that Xjv(O) = Xh{Q)/2{E - V{q{Q))). The 

constant factor 2{E — V{q{Q))) causes no serious problem, since we are interested in orthogonal relations between 
initial vectors. Moreover, it is an easy matter to see that X^+ii^) = grad-fir(O) = grad_ff(0) on account of the 
assumption gradl^(O) = at the initial point, where we note that gradX and gradff are taken with respect to 
metrics, gj and g^, respectively. To verify that the other initial vectors, X ^(0), a = 1, ■ • • , iV — 1, iV + 2, • • • , 2iV, are 
orthogonal to Xh{Q) and to gradi?(0), we use the following four facts: (i) Xi(0), • • ■ ,Xn-i{Q) G -Epjo) C\Hp(^Q), (u) 
Xjv+2(0), • • • , X27v(0) G i^p(o) n Vp(o), (iii) Hp(j^-^ and Vp(o) are orthogonal with respect to the Euclidean metric, as 
is seen from (I38|) . (iv) restricted to the subspaces i?p(o) and Vp(o), the Jacobi metric and the Euclidean metric are 
conformal to each other; 

g,,|p(0)(Xi,X2) = 2{E-V)gE\pio-,{Xi,X2), Xi,X2 G i/p(o), 

gj\p^o){X,,X2)=gE\pio)iXi,X2)/2{E-V), Xi,X2 G^Pfo). 

It then turns out from (i) and (iv) that Xi(0), • • ■ ,Xn-i{0) are also orthogonal to Xn{0) with respect to g^lpio), 
and further from (ii) and (iii) that they are also orthogonal to Xn+i{0) with respect to gE\p{o) ■ A similar statement 
for X]y+2{0), ■ ■ ■ ,X2n{0) holds true. 



D. A model system 



The model system we are to consider here is a natural Hamiltonian system with 3 degrees of freedom which has 
interactions of Henon-Heiles type. 



' 1 



i=l 



-p^ + VHHiq\q'+'] 



VHH{x,y) = x^y - iy^. 



(39) 



where = q^ . Hamiltonian vector fields both in the geometric method and in the usual method, denoted by Xk 
and Xh, respectively, arc given by ^ with g*^ = 5'^ /2{E - V) and V = X^Li 9*"^^)- 

Hamiltonian fiows of Xh for the Hamiltonian (|39|l are numerically calculated by the use of the 4-th order symplectic 
integrator which is a numerical integration method on the basis of discrete time evolution with each step 

an explicit symplectic mapping. Initial conditions for Hamilton's equations of motion are set as q^{Q) = and 
k(0) = Q:7i (j = 1, 2, 3), where 7i's are random values obtained from the uniform distribution function on the interval 
[0,1], and the constant a is determined so as to satisfy the energy condition = E. For the initial 

values ^'(O) = 0, we verify easily that the condition grady(O) = is satisfied, which was assumed in the previous 
subsection IV CI To integrate the linearized Hamilton's equations of motion ((SJ, we take an alternative method, that 
is, we choose to linearize, along a certain Hamiltonian flow, the sequence of symplectic mappings already obtained on 
the symplectic integrator algorithm. To our knowledge of explicit symplectic integrators, the symplectic integrator 
used here and another symplectic algorithm proposed in ^3] are set up on the assumption that Hamiltonians are of 
the form H{q,p) = T{p) + V{q), so that those algorithms are not applicable to the numerical integration of Hamilton's 
equations of motion with Hamiltonians of the form K{q,p) = X]i=i pf/^i^ ~ This means that we have to take 

another algorithm to integrate Hamilton's equations of motion in our geometric method for Lyapunov analysis. What 
we use in this article is an implicit but symplectic 6th-order Runge-Kutta method (Kuntzmann & Butcher method 
[Tsj l. see Appendix El- However, we have to note here that we do not need to apply that Runge-Kutta method to 
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integrate numerically Hamilton's equations of motion for K, since the solutions to Ea. (|21|l coincide with Hamiltonian 
flows already obtained by the explicit symplectic integrator up to the parameter change. We apply the implicit Runge- 
Kutta method to the numerical integration of the lifted Jacobi equations %V2\ . the linearized Hamilton's equations of 
motion for K. The implicit Runge-Kutta method, however, requires an additional process of numerical computation. 
In fact, we need to calculate the inverse of a 6A^ x QN matrix at each step of the integration, where N denotes the 
degrees of freedom. For this reason, the CPU-time we have needed to integrate the lifted Jacobi equations by the 
implicit Runge-Kutta algorithm is about 26 times as long as the CPU-time we have needed to integrate the linearized 
Hamilton's equations of motion for H by the explicit symplectic integrator. We have set the unit time slice as wide 
as /i = 2.5 X 10^^ both for the explicit symplectic integrator and the Runge-Kutta algorithm. 



E. Results of numerical calculations 




-0.01 ^ ' ' ' ' 1 

1000 2000 3000 4000 5000 

t 



FIG. 1: Convergence of Lyapunov exponents with E = 0.04. Curves represent graphs of Ka and Aa (a — 1,2,3), where 
and Aa, functions in the time parameter, are obtained by the geometric method and by the usual method, respectively. 



Figure n shows that Lyapunov exponents calculated in both methods have indeed definite values for E = 0.04, 
where A^'s and Ai^'''s are defined, respectively, to be 

l.WMm . 1, \\Va{t)\\ 

which are supposed to be convergent to Lyapunov exponents; limj^oo Ao(i) = Aa and limt^oo Ai^^(i) = Here, 
the quantities with the superscript (g) are those used in the geometric method. However, to compare the numerical 
results, we have made K^a\s) into a function of t by means of the parameter change. It is to be noted here that 

K'f \s) always vanishes on account of the fact that = \\Xk\\ = 2is:=constant. For E = 0.01, 0.02 and 0.03, 

we have obtained also definite Lyapunov exponents, which are shown in Fig|21 along with the dependence on energy. 
Figure 13 also shows that the Lyapunov exponents, Aa and Ai^\ calculated in both methods coincide with each other, 
which means that the Lyapunov exponents are obtained independently of the choice of methods, geometric or usual. 

We remark here that if one uses the Jacobi equations a second-order differential equation, to calculate the 
exponential growth rates of trajectories, one may obtain the same value as that obtained in the usual method. For 
example, for the Fermi-Pasta-Ulam (3 model, the largest Lyapunov exponent is calculated by using a 2iV-dimensional 
vector (X*, AX"^ /At) T^, where {X^{t)) is a solution to the Jacobi equations © and the Euclidean metric is used for 
the 2A^-dimensional vector. According to ,19j . the resultant value of the exponent coincides with the largest Lyapunov 
exponent obtained in the usual method. This might suggest that to calculate the largest Lyapunov exponent, one 
does not need to work with the cotangent spaces. However, the advantage of the geometric method developed in this 
article is that after the geometric method, we can obtain all the Lyapunov exponents along with the Lyapunov vectors 
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FIG. 2: Comparison of Lyapunov exponents obtained by both the geometric method and the usual method. By Aa and Xa, 
we denote Lyapunov exponents obtained by the geometric method and by the usual method, respectively. Numerical results 
obtained in both methods are in good agreement. 



among which two vectors associated with the vanishing Lyapunov exponents can be separated out from the others. 
This can be observed in Figs. 3 and 4. 

From FigsEland^ we will observe that the Lyapunov vectors calculated numerically in the geometric method satisfy 
the requirements stated in Sec and that the Lyapunov vectors calculated in the usual method have the property 
shown in Sec lVBI 

Figure 13 provides temporal evolutions of inner products between normalized Lyapunov vectors and the normalized 
tangent vector to a Hamiltonian flow. The inner products both in the geometric method and in the usual method are 
denoted by ii^^ and by ta, respectively: 



Figure 13 shows that all the Lyapunov vectors except for in the geometric method are orthogonal to Xk, and 

that the normalized V^^^ is equal to Xk- On the other hand, we observe also from FigOthat Vi and V2 in the 
usual method are not always orthogonal to Xh, and that the normalized V3 does not equal 1|, either. In 

particular, we remark that ^2 takes values around unity in opposition to our requirement. 

Figure ^ provides temporal evolutions of inner products between normalized Lyapunov vectors and the normalized 
gradient vector of the Hamiltonian, and the inner product both in the geometric method and in the usual method is 
denoted by the symbol ni®^ and Ua, respectively: 



All the Lyapunov vectors except for V/ are observed to be orthogonal to gradi^T, and Vf to be collinear to gradii' 
in the geometric method, as is expected. On the other hand, the Lyapunov vectors V^s, Vg in the usual method are 
not always orthogonal to gradi/, and V4 does not point to the direction of gradi/, either. In particular, 715 in the 
usual method is far from vanishing, taking values around minus unity. 

These observations agree to what we expect from the theory described in Sees II V 51 and IVBl We note in conclusion 
that tiny fluctuations around straight lines in Figs[3|and0I in particular, Figs|3Jb) and^e), seem to be numerical 
errors due to the factor 1/2 (i? — V) included in the metric g*-' , Christoffel symbol r*^, and the Riemann curvature 





tensor Rt-jki- 
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FIG. 3: Temporal evolutions of inner products between the normalized tangent vector of a Hamiltonian flow and normalized 
Lyapunov vectors. The energy is set aX E = 0.04. In (a),(b), and (c), straight lines are graphs of ti^-* against the time 
parameter in the geometric method, and broken curves are from the usual method, providing the graphs of ta- The 1st and 
2nd Lyapunov vectors V^-{'\ V'f^ are always orthogonal to the tangent direction to a Hamiltonian flow, Xk, in the geometric 
method, but Vi, 'V2 are not always orthogonal to Xh in the usual method. Moreover, the 3rd Lyapunov vector always points 
to the direction of Xk in the geometric method, but does not point to the direction of Xh in the usual method. In (d),(e), 
and (f), only straight lines are drawn, which are graphs from the both methods, but they coincide with each other. 



VI. CONCLUDING REMARKS 



In this paper, we have developed a new geometric method in Lyapunov analysis for natural Hamiltonian systems 
with N degrees of freedom, which is set up on the cotangent bundle of a Riemannian manifold endowed with the 

Jacobi metric. In contrast with our geometric method, the old or already-known geometric method is established on 
the Riemannian manifold with the Jacobi metric. According to that method, one brings Newton's equations of motion 
for a natural dynamical system into geodesic equations for the Jacobi metric and uses Jacobi equations, linearized 
geodesic equations, to analyze orbital instability of trajectories. However, the Jacobi equations are second-order 
differential equations, while Lyapunov exponents and vectors are defined through first-order differential equations. 
We then need first-order differential equation to apply Lyapunov analysis. According to our method, the Jacobi 
equations are lifted from Riemannian manifolds to their cotangent bundles to take the form of first-order differential 
equations. 

When the geometric method is applied, a question arises as to whether Lyapunov exponents remain unchanged in 
their values or not, in comparison with those obtained in the usual method. As we have already pointed out, the 
linearized equations in both methods are different from each other and can not be transformed to each other through 
the parameter transformation ds = 2{E — V{q))At, while the equations of motion in both methods are transformed 
to each other through the same parameter transformation. However, the numerical computation has shown that 
the values of Lyapunov exponents coincide with each other, independently of the choice of methods applied, as far 
as the model system with 3 degrees of freedom is taken. We guess that the Lyapunov exponents are long-term 
averaged values, so that they are independent of the choice of Lyapunov vectors along trajectories, while Lyapunov 
vectors depend on the choice of methods. As for the parameters of trajectories in the both method, we assume that 
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FIG. 4: Temporal evolutions of inner products between the normalized gradient vector of the Hamiltonian function and the 
normalized Lyapunov vectors. The energy is set at _E = 0.04. In (d), (e), and (f), straight lines are graphs of gainst the 

time parameter in the geometric method, and broken curves are from the usual method, providing the graph of n^. The 4th 
Lyapunov vector V'f'' always points to the gradient direction of the Hamiltonian function K, but V4 does not always point 
to gradiif in the usual method. Moreover, the 5th and 6th Lyapunov vectors Vg*', Vg^' are always orthogonal to gradif in 
the geometric method, and VsjVa are not so to grad^f in the usual method. In (a),(b), and (c), straight lines from the two 
methods are drawn, but each of them looks like a single line because of coincidence. 



the change of parameters must be subject to the condition < ds/dt < 00 along trajectories. On this account, 
we expect that Lyapunov exponents are independent of the choice of methods for calculation. We will find indeed 
the coincidence of Lyapunov exponents in the both methods from numerical computations for other model systems. 
Further, observations made from the Lyapunov exponents are expected to be independent of the choice of methods. 
For instance, a characteristics of the graph of Lyapunov spectra against i/N, i = 1, • ■ • ,N [2(1 l2ll |. which are 
observed in the usual method for a wide class of Hamiltonian systems having nearest neighbor interactions, will be 
found, in the geometric method as well, to be the same as that observed already in the usual method. We wish our 
geometric method may afford a fresh insight into the observation through Lyapunov vectors. 

In our geometric method developed in this article, we can choose Lyapunov vectors so as to satisfy the following 
requirements: (i) Lyapunov vectors except for A^-th and (N+ l)-th vectors are always orthogonal to both the tangent 
direction to a trajectory and the gradient direction of the Hamiltonian function, (ii) A'^-th Lyapunov vector points to 
the tangent direction of the trajectory, Xk-, and (iii) (N + l)-th Lyapunov vector points to the gradient direction of the 
Hamiltonian function, gradi^T. Along with such Lyapunov vectors, we can analyze orbital instability of Hamiltonian 
flows in phase spaces without influence of the two marginal directions pointed by Xk and gradif which have vanishing 
Lyapunov exponents, Xn = Xn+i = 0. Moreover, the iV-th and the (TV + l)-th local Lyapunov exponents, which are 
averages of exponential growth rate in finite time, vanish on any time interval. The local Lyapunov exponents in the 
usual method are used, for instance, to distinguish nearly integrable systems from the others ^22]. 

In this article, we have considered the Hamiltonian function of the form H{q,p) = ^ J2ij S^''PiPj+V{Q) and developed 
the geometric method in Lyapunov analysis. However, the geometric method can be established for Hamiltonian 
functions of the form H{q,p) — i J2ij ^^''{<l)PiPj + ^(9)7 where (a*^ (q)) is the inverse of a metric tensor {aij{q)). In 
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this case, the Jacobi metric is defined to be gij{q) = 2{E — V{q))aij{q), and geodesic equations for this metric 

dV _j_p» dq^ dq^ _ ^ 
ds2 ^'^ ds ds 

prove to be equivalent to Newton's equations of motion 




dq^dq^^ _ _ dV 
jk\ dt dt ~ " dqi 



with the total energy fixed at E, where r*^, and < > are the ChristofFel symbols formed from the metric gij and 

\jk J 

aij, respectively, and s is the length parameter for the Jacobi metric g^, which is related to the parameter t by 
ds 

— ~ 2(E — V(q)). The geometric method we have developed in Lyapunov analysis of linearized Hamilton's equations 
dt 

of motion on the cotangent bundle is independent of the choice of the Riemannian metric chosen, so that the theorem 
stated in Sec lIVBl holds also true in this case. Hence, we may find Lyapunov vectors which satisfy the requirements 
mentioned frequently. 
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APPENDIX A: GEOMETRY OF COTANGENT BUNDLES 



Vector fields and Levi-Civita connection on a Riemannian manifold M are lifted to the cotangent bundle T*M, and 
thereby the relation between geodesies on M and geodesic flows on T*M will be made clear in geometric fashion. 
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1. Lift of vector fields on M 

The cotangent bundle T*M is endowed with the standard one- form 6, which is expressed locahy as 6* = ptdx^ . 
Note that the 9 is defined globally on T*M . This can be seen from the coordinate transformation on the non-empty 
intersection (|19|l . The exterior derivative of 6*, w := d0, is the standard symplectic form on T*M . 

For vector fields on M, a way to lift them is not unique. A canonical way is given as follows: For Y e SG{M), the 
lifted vector field Y is defined through the conditions 

= Y, CyO = 0, (Al) 

where tt* is the differential of the canonical projection tt, and C denotes the Lie derivation. For Y = Y^di, a 
straightforward calculation shows that the Y is put in the form 

d dY^ d 

^^^^ 

Furthermore, owing to Cartan's formula, C-yO = d{L{Y)9) + L{Y)d9, along with t{Y)9 — 9{Y), the latter of the 
conditions IjAip implies that —d{9{Y)) — l{Y)uj, which then shows that the Y becomes the Hamiltonian vector field 
associated with F := 9{Y) ~ PiY"^ . Thus one has 

r = = (A3) 
dpi dx^ dxi dp^ 

With respect to the adapted frame, the canonical lift Y takes the form 

Y = YW,-p^V,YWj, (A4) 

where 

dY^ 

In addition to the canonical lift, one can define another lift; for a vector field Y = Y'^di on M, the horizontal lift 
of Y is given on T*M by 

Y^ = YWi. (A5) 

From the transformation rule H25|) . the horizontal lift is shown to be defined independently of the choice of adapted 
frames. 

A curve x{t) in M is also lifted horizontally; a curve x^{t) in T*M is called a horizontal lift of x(t), if 7r{x'^{t)) — x{t) 
and if the tangent vector to x'* (t) is horizontal. To give an example of horizontal lifts of curves, we consider a geodesic 
x{s) with s the arc length parameter. Let ^(s) denote its tangent vector and let Pi{s) = gij^-'i-s)- Then a curve 
{x{s),p{s)) in the cotangent bundle T*M is shown to be a horizontal curve. In fact, differentiation of {x{s),p{s)) with 
respect to s along with the geodesic equation for x{s) provides 

^(x^(s),p,(s)) = iCr'^.Pk^^) = (A6) 

as is wanted. From Ea. (|A6|l along with ^*(s) = g'^^pe{s), we observe that the curve {x{s),p{s)) is a geodesic flow, an 
integral curve of X^: (see fI7\i ). 

Now we assume that ^ is a tangent vector field to a congruence of geodesies in M. According to (|A5|I . we can define 

~h 

the horizontal lift ^ on T*M. With the restriction pi = gij^^ {xj imposed, the Hamiltonian vector field Xk becomes 

~h 

equal to the horizontal lift ^ . Hence, a congruence of geodesies in M is lifted to a family of geodesic flows in T*M 

~h 

along with ^ — Xk ■ 

We proceed to discuss lifts of geodesic deviations. Let Y{s) be a vector field defined along the geodesic x{s). We 
define a vector field X{s) along a geodesic flow {x{s),p{s)) with Pi{s) = gij£,^{s), by 

X = YW,+g,j{\7^Y)Wj. (A7) 

We note here that the X{s) is defined independently of the choice of adapted frames. If X{s) satisfies the lifted 
Jacobi equation H32|) . then the Y{s) should be a Jacobi field. Conversely, for a given Jacobi field Y{s) defined along 
a geodesic x(s), we can form a lifted vector field X{s) according to ljA7p . which is defined along a geodesic fiow 
P(s) = {x{s),pis)) with pi{s) = grj^^s). Then X{s) solves Eq.JSS- 
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2. Killing vector fields 

We now wish to investigate the relation between the canonical lift l|A3|l and the lift l|A7|l . where ^ is viewed 
as the tangent vector field to a congruence of geodesies in M. To this end, we first consider symmetry of our 
Hamiltonian system with the Hamiltonian function K. We assume here that for a vector field 1^ on M the function 
F = 0{Y) = y^pi is a constant of motion; Xk{F) = —{K,F} = 0, where {•, •} denotes Poisson bracket. Then one 
obtains [Xk, Xp] = —X[x,f} ~ 0. This implies that Y = Xp satisfies the linearized equation \3'2li along any geodesic 
flow. On the other hand, the condition Xk{F) = holds, if and only if y is a Killing vector field, an infinitesimal 
isometry, i.e., CyQ = 0, as is easily seen. It is well known that every Killing vector field satisfies the Jacobi equation 
along any geodesic. 

Now we assume further that we are given the tangent vector field ^ to a congruence of geodesies in M. If restricted 
on a subspace L determined by pi — gij£_^ in T*M, the canonical lift Y of a Killing vector field Y is expressed as 

Y\l = Y'D, ~ g.k^'^V^YWj = Y'D, + g,jiV^Y)Wj, (A8) 
where use has been made of the formula that 

g'^V.^Y'' + g'^ViY^ = 0, (A9) 

which is a necessary and sufficient condition for I' to be a Killing vector field. Thus we have found that if V is a 
Killing vector on M, and if the canonical lift Y is restricted to L determined by pi = gij^^ , then Y\L is equal to the 
lift HA7|) with ^ the tangent vector field to a congruence of geodesies. 

3. Levi-Civita connection of T*M 

The Levi-Civita connection V is defined on the cotangent bundle T*AI through the Sasaki metric g. We denote 
the Christoffel symbols for this connection by Tg^; 

V9,dc=ticdA, 

where Roman capital indices run from 1 to 2m and Oa are the standard frame; 

The Christoffel symbols are given, as usual, by 

^Bc = ^9^" {dsgcD + dcgoB - dogBc) , 

where gAB are components of g; gAB — g{dA, ds)- We denote the coefficients of the connection V with respect to the 
adapted frame by F^^; 

where Greek indices also run from 1 to 2m, but indicating that they are indices for the adapted frame. 
Let the functions f^^^" be defined by 

Then the torsion-free condition for V is put in the form 

A straightforward calculation yields f^^-y" as follows: 

[D,,D,] = PiR,,/D^, [A, Dj] = -ri,Dj:, 
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Wc are to write out F^^ in terms of ga/3 and f^^-y", where ga/s = g{Da, Dp), the components of g with respect to 

the adapted frame. The covariant derivative of the metric g must vanish for all vector fields X on T*M; Vxg = 0, 
so that one has 

D^g^s - r^^5£5 - r^^^T-e = 0. 

Further calculation provides 

Dpg^s + D^gsp - Dsgp^ = (f^^ + f g,s + (f ^5 - 9e^ + {^6 ' ^h) 3ep 



which results in 



where 



F^^ = \r' {Dffgjs + D,gs0 - Dsgp^) + \ + + 



fS-y — 9 ''''6/3 9ei- 

A straightforward calculation shows that the components F^^ are given by 



f}. = r}„ f;^ = -ii?fV, ri, = -li?;^V, r^ = o, 



^k = lRjkM r5-, = -F^, F|, = 0, F^ = 0. 
Covariant derivatives of vector fields are then expressed, in terms of these coefficients, as 



(AlO) 



Da, (All) 



where (Xf ) and (X^) are components of Xi and X2 with respect to Da, respectively. In particular, the covariant 

derivative of X = X^Di + X^Dj with respect to the horizontal lift ^ = £,^{s)Di along a geodesic flow as a horizontal 
lift of a geodesic takes the form 

~ H 1 — 

(V ,x)^ = 'i^+Ti-ex^ - ^Ri'%ex^, 

^ (A12) 
~ - H - 1 



~h 

If X = ^ , these equations give rise to 



v~^Hi = 0, 



which implies that the horizontal lift {x{s),p{s)) of a geodesic x{s) on Af is also a geodesic on T*M with respect 
to the lifted metric g. We note here that the arc length parameter a with respect to g reduces to the arc length 
parameter s, if the curve is horizontal. 

APPENDIX B: SYMPLECTIC IMPLICIT RUNGE-KUTTA METHOD 

Suppose we are given a dynamical system in 

d T 

-^^{t) = f{x,t). (Bl) 
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TABLE II: Kuntzmann & Butcher method, order 6. The upper right block is the matrix {aij), the lower raw is the vector [bi) 
and the left column is the vector (ci). 



1 

2 


~ nr 




5 

36 


2 
9 


^/T5 
~ T5~ 


5 

36 


~ "30" 




1 

2 


5 

36 


^/T5 




2 
9 


5 

36 


~ "24~ 


1 
2 


VT5 

10 


5 

36 


\/T5 
30 


2 
9 


\/T5 
15 




5 

36 






5 

18 




4 
9 




5 

18 



Numerical integration of this equation is performed through discretizing it with time shoe h. The s-stage Runge-Kutta 
method for integration is given by 

s 

x' = X + h biki 
1=1 

s 

h = fix + h'^a.ijkj^t + ah), i = 1, • • • , s, 
i=i 

where {x,t) goes to {x\t + h) after one step, and aij, hi and Ci are real constants with X]i=i = 1- Note that 
the second of the above equations defines imphcitly ki. The 3-stage Runge-Kutta method, namely the 6-th order 
Kuntzmann & Butcher method, is defined as in Table ITTl 



